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ABSTRACT 


Several  methods  for  reconstructing  the  resistivity  profile  of  a  layered, 
laterally  homogeneous  earth  from  direct  current  measurements  are  described. 

These  methods  recover  the  resistivity  of  the  earth  layer  by  layer  in  a 
recursive  way,  and  require  a  very  small  amount  of  computational  effort.  They 
are  obtained  by  transforming  the  inverse  resistivity  problem  into  an  equivalent 
inverse  scattering  problem,  and  by  applying  efficient  signal  processing  algorithms 
such  as  the  Schur,  fast  Cholesky  or  Levinson  recursions  to  the  transformed 
problem.  These  algorithms  operate  on  a  layer  stripping  or  layer  accumulation 
principle,  and  are  shown  to  be  related  to  previous  reconstruction  techniques 
of  Pekeris,  Koefoed,  Kunetz  and  Rocroi,  and  others. 
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1 .  Introduction 

The  problem  of  reconstructing  the  resistivity  of  a  layered  earth 
model  from  direct  current  measurements  has  been  the  object  of  sustained 
interest  over  the  years.  In  this  problem,  some  direct  current  is  injected 
inside  the  earth  by  two  current  electrodes,  and  two  voltage  electrodes  are 
used  to  measure  potential  variations  on  the  surface  of  the  earth.  The  goal 
is  to  reconstruct  the  resistivity  profile  as  a  function  of  depth  from  the 
potential  measurements  on  the  surface  of  the  earth.  The  existence  of  a 
solution  for  this  problem  was  established  by  Slichter  [1]  and  Langer  [2] , 
whose  solution  was  however  impractical  from  a  computational  point  of  view. 

In  1940,  Pekeris  [3]  obtained  some  recursions  for  reconstructing  the  earth 
layer  by  layer.  This  reconstruction  method  was  subsequently  refined  and 
developed  more  fully  by  Koefoed  [4] ,  [5] .  More  recently,  Coen  and  Yu  [6] 
used  a  transformation  procedure  of  Weidelt  17]  to  formulate  the  inverse 
resistivity  problem  of  the  earth  in  such  a  way  that  the  Gelfand-Levitan 
method  of  inverse  scattering  theory  could  be  applied  to  this  problem.  Another 
layer  by  layer  reconstruction  method  was  also  proposed  by  Kunetz  and  Rocroi  [8] , 
and  it  will  be  shown  below  that  their  algorithm  can  be  identified  with  the 
Levinson  recursions  of  linear  prediction  [9] ,  which  in  fact  arise  in  a  large 
number  of  signal  processing  situations. 

The  objective  of  this  paper  is  to  give  a  unified  account  of  layer  by 
layer  reconstruction  techniques  for  the  earth  resistivity,  which  includes 
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previous  reconstruction  methods ,  and  new  ones  as  well .  The  common 
framework  that  will  be  used  for  this  presentation  is  that  of  inverse  scattering, 
and  in  this  context,  we  will  show  that  layerwise  reconstruction  procedures 
can  be  viewed  as  differential  inverse  scattering  methods  of  the  type  discussed 
by  Bube  and  Burridge  [10] ,  and  Bruckstein,  Levy  and  Kailath  [11]  (see  also 
[12]  -  [14]).  Since  the  equation  describing  the  potential  of  the  earth  is 
elliptic,  the  inverse  resistivity  problem  of  the  earth  is  not  an  inverse 
scattering  problem.  However,  by  writing  the  equations  for  the  earth's  potential 
as  two  coupled  first-order  equations,  we  will  be  able  to  introduce  a  matrix 
whose  elements  can  be  viewed  as  obtained  by  analytic  continuation  of  the  elements 
of  the  scattering  matrix  associated  to  a  true  scattering  system.  In  this 
context,  we  show  that  the  transformation  of  Weidelt  [7]  and  Coen  and  Yu  [6] 
has  for  effect  to  map  solutions  of  the  potential  equation  into  solutions  of 
a  wave  equation  whose  scattering  matrix  is  the  one  mentioned  above. 

The  advantage  of  formulating  the  inverse  resistivity  problem  as  an 
inverse  scattering  problem  is  that  the  relation  between  layerwise  reconstruction 
methods  (which  are  also  called  layer  stripping  techniques)  for  this  class 
of  problems,  and  efficient  signal  processing  algorithms  such  as  the  Schur, 
fast  Cholesky  and  Levinson  recursions  has  been  the  object  of  close  scrutiny 
[10] ,  [11] ,  [15] .  We  will  show  for  example  that  the  recursions  obtained  by 

Pekeris  [3]  and  Koefoed  [4] ,  [5]  for  reconstructing  the  earth  resistivity 

are  just  a  modification  of  the  Schur  algorithm  [16] ,  [17]  which  is  now  widely 
used  in  linear  estimation  theory,  or  network  synthesis  [18] .  The  continuous 
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parameter  version  of  this  algorithm,  which  takes  the  form  of  a  Riccati 
equation,  will  also  be  related  to  the  work  of  Langer  [2]  and  Slichter  [1] . 

In  addition,  it  will  be  shown  how  efficient  algorithms  such  as  the  fast 
Cholesky  or  Levinson  recursions  can  be  used  to  recover  the  earth's  resistivity. 
The  solution  based  on  the  Levinson  recursions  turns  out  to  be  identical  to 
the  method  proposed  by  Kunetz  and  Rocroi  [8] . 

This  paper  is  organized  as  follows.  In  Section  2,  the  inverse  resistivity 
problem  is  described  and  its  relation  with  an  equivalent  inverse  scattering 
problem  is  examined.  Section  3  describes  the  solution  of  the  inverse 
resistivity  problem  via  the  Schur  algorithm  and  its  relation  to  the  work  of 
Pekeris,  Koefoed  and  Langer.  In  Section  4,  it  is  shown  that  after  applying 
a  transformation  similar  to  that  of  Weidelt  [7] ,  the  fast  Cholesky  and  Levinson 
recursions  can  be  used  to  recover  the  resistivity  of  the  earth.  The  relation 
of  these  methods  with  those  of  Coen  and  Yu,  and  Kunetz  and  Rocroi  is  also 
discussed.  In  Section  5  we  describe  how  the  given  data,  which  is  usually 
the  apparent  resistivity  for  the  Schlumberger  electrode  configuration  [5] , 
can  be  used  to  compute  the  functions  used  to  perform  the  inversion  with  the 
Schur  algorithm  or  the  fast  Cholesky  and  Levinson  recursions,  which  are 
respectively  Slichter' s  kernel  function  [1],  [5]  and  a  certain  fictitious 
current  source  profile  obtained  by  Maxwell's  method  of  images  ([8],  [5], 

Chapter  10) .  Finally,  Section  6  contains  some  conclusions  and  some  suggestions 


for  further  research. 
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2.  Problem  Formulation 

The  inverse  resistivity  problem  of  the  earth  is  formulated  here  as 
in  [1]  -  [6].  It  is  assumed  that  the  earth  conductivity  0( z)  varies  with 
depth  only,  and  in  a  first  stage  we  consider  the  idealized  problem  where 
some  direct  current  flows  inside  the  earth  through  a  single  electrode,  and 
where  the  potential  <j>(0,r),  where  r  denotes  the  radial  distance  to  the 
current  electrode,  is  measured  on  the  surface  of  the  earth.  The  objective 
is  to  reconstruct  0( z)  from  <f>(0,r).  In  a  second  stage,  we  will  examine 
the  more  realistic  situation  where  the  Schlumberger  electrode  configuration 
is  used  to  measure  the  apparent  resistivity 

pa(r)  =  “  y-  r2  4>(0,r)  ,  (2.1) 

where  <j)(0,r)  is  the  potential  obtained  for  a  single  current  electrode,  and 
I  is  the  current  supplied  by  the  source.  The  objective  in  this  case  will 
be  to  reconstruct  0(z)  from  p  (r), . 

cl 

When  a  single  current  electrode  is  used,  by  selecting  this  electrode 
as  the  origin  of  cylindrical  coordinates  (z,r,0),  the  current  equations  are 

j(z,r)  =  0(z)  V<j>(z,r)  (2.2) 

V*j(z,r)  =  0  (2.3) 

where  <j>(z,r)  is  the  potential,  and 


jZ (z,r)" 

j(z,r)  = 

-j  (z,r)_ 


(2.4) 
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is  the  decomposition  of  the  current  density  in  its  vertical  and  lateral 
components.  The  boundary  conditions  for  these  equations  are 


.2  .  I  6(r) 

3  (0,r)  =  -  — 


21T  r 


(2.5) 


where  I  is  the  total  current  supplied,  and  $(z,r)  0  and  j(z,r)  -*■  0  as 

2  2  1/2 

R  =  (z  +  r  )  -*•  00 .  By  combining  C2.2)  and  (2.3)  we  obtain  the 


potential  equation 


■q 

Acf>(z,r)  +  —  £n  a(z)  v—  4>(z,r)  =  0 
dz  a  z 


(2.6) 


3^  1  3  32 

where  A  =  Laplacian  =  — —  +  —  q—  +  — —  ,  which  is  the  equation  usually 

«  2  r  dr  ^  2 

dz  dr 

used  to  analyze  the  inverse  resistivity  problem.  However,  instead  of 
focusing  our  attention  on  equation  (2.6),  we  will  use  here  the  first-order 
equations  (2.2)  -  (2.3)  to  formulate  the  inverse  resistivity  problem. 

Let 

r 

(2.7) 


n 


H  [f(r)]  =  f  (x)  J  (Ar)rdr 


0 


be  the  Hankel  transform  of  order  n  of  a  function  f ( • ) ,  and  denote 


<}>(z,A)  =  HQ  [(f>(z,r)  ] 
j2  (z, A)  =  HQ [ jZ (z,r) ] 
jr(z,A)  =  H],[jr(z,r)  ]  . 


(2.8a) 

(2.8b) 

(2.8c) 


By  using  the  property 


H  [—  f  (r)  +  —  f  (r)  ]  =  Ah  [f  (r)  ] 
Or  dr  1 


(2.9) 
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Hi[^  f(r)]  =  _  ^Vf  {r)] 

of  Hankel  transforms,  the  equations  (2.2)  -  (2.3)  become 
jr(z,A)  =  -  Aff(z)  <j>(z,A) 
jZ(z,A)  =  a(z)  ^  $(z,A) 

Ajr(z,A)  +  jZ(z,A)  =  0 


<?r 


Eliminating  j  ,  and  denoting 


p,  -  -v ,  A 
ip(z,A)  = 


1  ^z 

^  j  (z,A)  = 


0  ( z )  d 


A 


dz 


<J>(z,A) 


this  gives 


(2.10) 

(2.11a) 

(2.11b) 

(2.12) 

(2.13) 


_d_ 

dz 


<f>(z,A) 

Jp( z ,  A)_ 

- 

Aa  1(z)" 


<P(z,A) 


.ip(z,A). 


which  is  the  analog  of  the  telegrapher ' s  equation 


d 

v(z,Af  0  -jAz(z)- 

v(z, A) 

dz 

/v  -  -1 

_i(z,A)J  (z)  0 

_  i^ZjA), 

(2.14) 


(2.15) 


satisfied  by  the  voltage  and  current  along  a  nonuniform  transmission  line, 
which  was  the  starting  point  of  the  inverse  scattering  problem  considered 
by  Bube  and  Burridge  [10],  and  Bruckstein,  Levy  and  Kailath  [11].  Note  however 
that  there  is  an  important  difference  between  (2.14)  and  (2.15):  A  is 
replaced  by  jA.  This  difference  arises  from  the  fact  that  the  equation 
satisfied  by  <f>  is  elliptic,  whereas  the  equation  satisfied  by  the  voltage 
along  a  nonuniform  transmission  line  is  hyperbolic. 
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This  difference  prevents  us  from  formulating  the  inverse  resistivity 
problem  as  an  inverse  scattering  problem,  since  inverse  scattering  theory 
applies  only  to  hyperbolic  operators.  However,  by  using  a  mapping 
technique  originally  introduced  by  Weidelt  [  7 ]  (see  also  Coen  and  Yu  [6]) , 
we  will  show  below  that  the  inverse  resistivity  problem  can  be  transformed 
into  an  equivalent  inverse  scattering  problem. 

A 

The  initial  conditions  for  the  differential  system  (2.14)  are  <j>(0,A), 
which  is  obtained  by  taking  the  Hankel  transform  of  the  observed  potential 
<p  (0  ,  r).  on  the  surface  of  the  earth,  and 

V(  0,AJ  =  .  (2.16) 

Following  1111  ,  we  introduce  the  normalized  variables 

M(z,A)  =  a1/2(z)$(z,X)  (2.17a) 

N(z,A)  =  0~1/2(z)H z,A)  (2.17b) 


so  that 


M(z,  A) 
N  Cz ,  A) 


a(z)  1^-  . 

\p(z,A) 


Then,  if  the  down  and  upgoing  waves  are  defined  as 


D  ( Z , A )  =  y  (M(z,A)  +  N(z,A)) 
D(z,A)  =  ^  (M(z, A)  -  N(z,A) ) , 


the  system  (2.14)  can  be  rewritten  as 


(2.18) 

(2.19a) 

(2.19b) 


D (z , A)  -  A  -  k(z)  d(z,A) 

—  =  (2.20) 

U(z,A)  -  k(z)  A  U(z,A)_ 
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where  k(z)  =  -  X-  Unff  (z)  is  the  reflectivity  function.  By  discretizing 
(2.20),  we  obtain  the  elementary  filter  sections  described  in  Fig.  1.  These 
sections  show  that  the  waves  D(.z,A)  and  U(z,A)  propagate  in  opposite  down¬ 
ward  and  upward  directions,  and  for  a  layer  of  thickness  A  at  depth  z, 

D(z,A)  and  U(.z,A)  are  attenuated  by  a  factor  exp  (-AA)  and  are  partially 
reflected  in  the  proportion  k(z)A. 

We  have  that  k(z).  =  0  for  z<  0,  and  we  assume  that  k(*)  is  summable 
and  has  compact  support,  so  that  there  exists  L  >  0  such  that  k(z)  =  0 
for  z  >  L.  In  this  case  the  two-component  system  (2.20)  can  be  viewed  as 
perturbed  form  of  the  free  system 


_d_ 

'  Dq(z,A) 

A 

0 

Dq(z,A) 

dz 

VZ,A). 

0 

A 

_U0(z,A)_ 

(2.21) 


where  the  perturbation  k(')  is  small,  so  that  for  z  <  0  and  z  >  L,  the 
solutions  of  (.2.20)  are  identical  to  those  of  (2.21),  i.e. 


D  Cz ,  A)  =  D  (A) e  Az  (2.22) 

Xj 

U(z, A)  =  U  (A)e^Z 

L 


for  z  <  0,  and 

D (z, A)  =  D  (A)e~X 

K 

U(z,A)  =  U  (A)eAz 

K 


(2.23) 


for  z  >  L.  By  linearity,  we  have 
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r  d  (X)  i 

rD  (A)-i 

R 

=  s  (A) 

L 

u  (A) 

-L  “ 

.VA). 

where 


S(A)  = 


T  (A) 

r  (A)i 

L 

R  1 

R  (A) 

Xj 

te(A). 

(2.24) 


(2.25) 


is  the  scattering  matrix  of  (.2.20).  As  shown  in  Fig.  2,  this  matrix  relates 
the  incoming  and  outgoing  waves-  for  the  aggregate  medium  obtained  by 
composing  the  elementary  layers  described  in  Fig.  1. 

The  motivation  for  calling  S  (X)  a  scattering  matrix  is  that  by  setting 
X  =  jk  with  k  real  in  (2.20),  we  obtain  a  two-component  scattering  system 
of  the  type  discussed  in  119],  (20],  Ill],  [15],  The  scattering  matrix 
of  such  a  system  is  S(.jk).  It  has  the  property  that 


SH  ( jk)  S  (jk)  =  I 


(2.26) 


for  k  real,  where  the  superscript  H  denotes  the  Hermitian  transpose,  and 
it  obeys  the  reciprocity  relation 

T  (X)  =  T  (A)  =  T (X)  (2.27) 

J_i  K 

for  all  X.  In  addition,  since  k(*)  is  summable,  the  two-component  wave 
system  (2.20)  with  X  =  jk  has  no  bound  states  ([19],  Chapter  1),  and  T(X) 
is  analytic  in  the  right  half-plane.  Similarly,  k(z)  =  0  for  z  <  0  implies 
that  R  (X)  is  analytic  in  the  right  half-plane,  and  the  assumption  that  k(.) 
has  compact  support  implies  that  R  (X)  is  meromorphic  in  the  right  half-plane 


[21J  . 
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Thus,  R  (A)  and  T(A)  with  Re  X  >  0  can  be  viewed  as  obtained  by 
-Li  *“ 

analytic  continuation  of  R  (jk)  and  T(.jk) ,  whereas  the  assumption  that 

Li 

k(z)  has  compact  support  is  necessary  to  guarantee  the  existence  of  Rr(X) 

in  the  right  half-plane  (note  that  it  may  not  be  defined  at  some  points) . 

However,  since  our  analysis  below  will  focus  exclusively  on  R  (X) ,  this 

L 

assumption  may  easily  be  removed. 

For  the  problem  considered  here,  the  earth  is  probed  from  its  surface, 
so  that  U  (A)  =  0  in  (2.24)  and  the  reflection  coefficient  R_ (X)  can  be 

-K  _L 

expressed  as 


TT  (' ^  /s 

w  m  _  L  J  _  O'  (0)  $  (0,X)  -  Tfi( ,0,X) 

D  (A)  ^  ^ 

'  a  (,0)<j)  (o,A)  +  ip(.OrA) 


(2.28) 


where  <})(.0,X)  is  the  transform  of  the  observed  potential,  and  ijj(0,X)  is  given 
by  (.2.16).  In  the  resistivity  prospecting  literature  R  (X)  is  known  as  the 

Li 

modified  kernel  function  ([5],  p.  202).  It  is  related  to  Slichter's  kernel 
function  K(A)  by 


R  m  -  K(X).  ~  1 

VA)  K(X)  +  1  ' 


(2.29) 


where  K(X)  is  the  normalized  impedance  of  the  resistive  medium  extending 
over  IO,00),  i.e. 


KM  =|^4f  =0(0,  4^ 

N<0M  ifto.X) 


(2.30) 
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Since  K(A)  and  R  (A)  are  entirely  specified  by  the  given  data,  the 
L 

inverse  resistivity  problem  can  be  posed  as  follows:  given  K(A)  or  R  (A) 

for  A  real  and  positive,  we  want  to  reconstruct  k(z)  and  O {z) .  In  theory, 

this  can  be  done  by  using  the  fact  that  R  (A)  is  analytic  in  the  right 

Xj 

half-plane  to  obtain  its  value  on  the  imaginary  axis,  and  then  by  using 
any  of  the  inverse  scattering  techniques  described  in  [10],  [11]  to  recover 
k(z)  from  R  (jk)  .  However,  this  basic  scheme  can  be  implemented  in  a 
variety  of  ways,  which  we  will  now  discuss  and  compare. 
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3*  The  Schur  Algorithm  and  Related  Methods 


The  reconstruction  techniques  which  will  be  examined  in  this  section 
can  all  be  viewed  as  variants  of  an  algorithm  introduced  by  Schur  [16]  in 
1917  (see  also  Akhiezer  117] )  to  test  the  boundedness  of  an  analytic  function 
inside  the  unit  circle,  and  which  solves  the  inverse  scattering  problem  for 
a  discrete  medium  [18]  .  This  algorithm  was  subsequently  extended  to  con¬ 
tinuous  two-component  wave  systems  in  Ill.]  ,  [15]  ,  and  we  will  now  adapt 
this  version  of  Schur' s  algorithm  to  the  system  (2,20). 

We  denote  by 


r  ,  i  ,  U(z,A) 

VZ'A)  "  SteAT 


(3.1) 


the  reflection  coefficient  associated  to  the  section  of  the  resistive 

medium  extending  over  [z,00)  .  By  using  the  differential  equation  (.2.20) 

for  D(z,A)  and  U(z,A),  we  find  that  R  ( z ,  A )  satisfies  the  Riccati  equation 

L 


Rl(z,A)  =  2AR^(z,A)  +  k  (z)  (R^(z,A)  -  1)  , 
with  the  initial  condition 


(3.2) 


Rt  (0, A)  =  R  (A)  (given)  .  (3.3) 

i-i  L 

The  equation  (3.2)  depends  on  k(z),  and  therefore  if  we  want  to  use  it 
to  reconstruct  k(-),  we  need  to  express  k(z)  as  a function  of  R  (z,A).  To 
do  so,  note  that  since  R  (z,A)  is  the  analytic  continuation  of  R  ( z , j k ) 

Xj  L 

with  k  real,  where  R  (z,jk)  -*■  0  as  k  -*■  00  [11],  it  can  be  expanded  as 

Jj 


CO 

=  E  r-= 

*  ... 


Rl(z, A) 


(3.4) 
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so  that  after  substitution  in  (3.2) ,  we  obtain 


k(z)  =  2r  (z)  =  lim  2  A  R  (z,A)  ,  (3.5) 

1  A-*»  L 

By  substituting  (3.5)  inside  (3.2),  R (z,A)  can  be  propagated  recursively, 

L 

and  in  the  process  we  reconstruct  k(.z)  for  all  z.  This  gives  also  CT(z) 
by  noting  that 

fz 

0 (z)  =  a (0)  exp  (-2  k(s)ds)  (3.6) 

■'0 

where  O (0)  is  the  conductivity  on  the  surface  of  the  earth. 

The  recursions  (3.2),  (3.5)  constitute  the  Schur  algorithm.  This 

reconstruction  procedure  can  be  viewed  as  the  continuous  parameter  version 

of  the  method  proposed  by  Pekeris  in  [3]  (see  also  [5] ,  Chapter  10) .  To 

see  this,  assume  that  the  conductivity  function  O(z)  is  approximated  by 

a  piecewise  constant  function,  so  that  the  earth  can  be  viewed  as  constituted 

of  N  homogeneous  layers  of  thickness  t^  and  conductivity  O ,  1  <  i  <  N. 

Then,  by  discretizing  (.3.2)  and  assuming  that  the  thickness  of  every  layer 

2 

is  sufficiently  small  so  that  terms  of  order  t^  can  be  neglected,  we  obtain 


2At .  (R.  (A)  -  k. ) 

_  .  x  i _  1 

Ri+1  A}  “  e  (i  -  k.R. (A)) 

i  x 


(3.7) 


with  R  (A)  =  R  (A) ,  where 
U  L 


R.  (A)  =  rt 

1  I_! 


Ui  v) 


(3.8) 


is  the  reflection  coefficient  obtained  by  stripping  away  the  i  first  layers 

of  the  resistive  medium,  and  by  assuming  that  the  current  electrode  is 

fcl"l 

located  on  top  of  the  i+1  layer.  From  (3.5),  we  find  also  that 


a . 
1 


-  a 


i+l 


a. 

i 


+  ai+i 


2At. 

1  R±a>  , 


k. 

l 


lim  e 
A-ko 


(3.9) 
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or  equivalently 


i  R.  (X)  =  In  k.  -  2At. 
n  1  lx 


(3.10) 


as  X  -*  00 ,  which  is  precisely  the  formula  used  by  Pekeris  to  reconstruct  t^ 
and  k^.  From  (3.10),  we  see  that  as  X  -*  00 ,  the  function  R^(A)  can  be 
approximated  by  a  line  whose  slope  is  -2t^  and  whose  intersection  with 
the  vertical  axis  is  £nk^.  Consequently,  by  combining  (3.7)  and  (3.10), 
R^(A),  t^  and  k^  can  be  computed  recursively  for  1  <_  i  <_  N,  and  (3.9)  can 
be  used  with  the  initial  condition  a.  =  0(0)  to  obtain  0.  for  all  i. 

0  l 

Instead  of  propagating  the  reflection  coefficient  R  (z,A)  ,  we  could 

±J 

choose  to  propagate  the  normalized  impedance,  i.e.  Slichter's  kernel 
M (z, A) 


K  ( z ,  A )  = 


N(z , A) 


(3.11) 


which  is  related  to  R  (z,A)  by  the  relation 

Ij 


r  (  , s  K(Z,A)  -  1 

Vz'A)  -  surn 


(3.12) 


By  noting  that 


d 

dz 


-k(z)  A  I  rM(z,A)‘ 

-  A  k (z)J  LN(z , A) . 
we  find  that  K(z,A)  satisfies  the  Riccati  equation 


"M  (z ,  A)" 

.N(z,A)_ 

(3.13) 


~~  K(z,A)  =  -  2k(z)K(z,A)  +  A(K2(z,A)  -  1) 
dz 


(3.14) 


with  the  initial  condition  K(0,A)  =  K(A)  (given) ,  which  was  first  derived 
by  Langer  [2].  By  expanding  K(z,A)  as 


K(z, A)  =  i|0  K  (z)A 


-l 


(3.15) 
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we  can  identify  KQ(z)  =  1,  and 


k(z)  =  K  (z)  =  lim  X(K(z,X)  -  1) 
A-*00 


(3.16) 


By  combining  (3.14)  and  (3.16),  we  can  therefore  propagate  K(z,X)  recursively 
and  in  the  process  reconstruct  k(z).  The  difference  between  this  recon¬ 
struction  method  and  that  of  Langer  12]  is  that  Langer  did  not  recognize 


that  the  inversion  could  be  performed  recursively.  Instead,  he  showed  that 
k(.0)  =  K^(.O)  and  that  all  derivatives  k  ^  (.0)  can  be  expressed  in  function 
of  K_.  (.0)  for  j  <_  i+1,  which  by  using  the  Taylor  series  expansion 


k(z>  -  Jo  fr 


(3.17) 


implies  that  k(z)  can  be  reconstructed  from 


K(X)  =  1  +  .  E,  K.  (0) X 

1=1  l 


(3.18) 


which  is  the  given  data. 

An  even  better  inversion  procedure  which  can  be  used  to  reconstruct 

a (z)  directly  (instead  of  k(z))  is  to  consider  the  unnormalized  impedance 

/\ 

Z  (z ,  X)  =  a_1(z)K(z,X)  =  (3.19) 

i/T(z,X) 

which  was  called  the  "resistivity  transform"  by  Koefoed  ([4],  [5],  Chapter  3) 
Then  Z(z,X)  satisfies  the  Riccati  equation 


~  z(z,X)  =  X(0(z)z2(z,X)  -  a  1(z))  ,  (3.20) 

dz 

and  by  using  the  expansion  (3.15)  with  the  observation  that  KQ (z)  =  1, 
we  find  that 

a  1(z)  =  lim  Z (z , X) 

X-x» 


(3.21) 
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so  that  0  (z)  can  be  recovered  by  propagating  (3.20)  and  (3.21)  recursively. 
This  reconstruction  procedure  is  the  continuous  analogue  of  the  recursions 

Z.  (A)  —  tanhCAt.) 

Zi+1(A)  =  1  -  a.  tanh  (At.)Z. (A)  (3.22) 

1  11 

with 

a"1  =  lim  Z.  (A)  (3.23) 

l  .  l 

A-**> 

which  were  obtained  by  Koefoed  [4] ,  {5]  to  reconstruct  a  discrete  resistive 
medium  constituted  of  N  horizontal  homogeneous  layers  of  thickness  t^  and 
conductivity  (7  ,  1  <_  i  <_  N,  where  Z^(A)  is  the  impedance  obtained  by 
removing  the  i  first  layers  of  the  medium.  A  straightforward  discretization 
of  (3.20)  -  (3.21)  can  in  fact  be  used  to  obtain  the  recursions  (3.22)  - 


(3.23)  . 
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4.  Fast  Cholesky  and  Levinson  Recursions 


The  inversion  procedures  that  were  described  above  use  either  equations 
(3.2),  (3.5)  to  reconstruct  k(z),  or  equations  (3.20)  -  (3.21)  to  reconstruct 
CT  (z) .  Since  these  methods  are  variants  of  the  Schur  algorithm  considered  in 
[11],  they  suffer  from  the  same  limitations.  The  most  significant  of  these 
is  that  we  need  to  take  the  limit  of  2Arl(z,A)  or  of  Z(.z,A)  as  A  -*■  00 ,  which 
is  not  a  very  reliable  numerical  operation.  To  eliminate  this  difficulty, 
we  will  now  show  that  the  problem  can  formulated  in  a  way  such  that 
efficient  signal  processing  algorithms  such  as  the  fast  Cholesky  or  Levinson 
recursions  [  9] ,  122]  -  1 23]  can  be  used. 

To  do  so,  we  will  use  the  method  of  Weidelt  [  7]  (see  also  [6])  to 
convert  the  inverse  resistivity  problem  into  an  equivalent  inverse 

^  v\ 

scattering  problem.  The  key  step  is  to  view  the  functions  A(f> (z , A)  ,  Ai|>(z,A), 

V  V 

Ad(z,A)  and  Au(z,A)  as  Laplace  transforms  of  some  functions  $(z,t),  ^(z,t), 

V  V 

D(z,t)  and  U(z,t) ,  so  that  if 


L[f  (t)  ] 


A 


•oo 

f(t)  exp(-At)dt 

0 


(4.1) 


denotes  the  Laplace  transform  of  a  function  f ( • ) ,  we  have 


A$(z,A) 

=  L[<f>(z,t)  ],  Ai|T(z,  A) 

V 

=  L[lf>(z,t)  1 

(4.2a) 

Ad(z,A) 

V 

V 

=  L[D(z,t)  ],  Au  (z,A) 

=  L  [U  (z ,  t)  ]  . 

(4.2b) 

By  multiplying  (2.14)  and  (2.20)  by  A,  and  taking  inverse  Laplace  transforms, 


we  obtain 
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_3_ 

3z 


<j)(z,t) 

V 

= 

_3_ 

St 


a”1  (z) 


-  v 

<j)(z,t) 

V 

J 

_^(.z,t)_ 

(4.3) 


and 


_S_ 

Sz 


.V 


D  (z ,  t) 

V, 

= 

U  (z ,  t) 

3t 

k  (z) 


-  k  (z) 

St 


D(z,t) 

V 

U(z ,  t) 

(4.4) 


which  are  respectively  the  telegrapher  and  two-component  wave  equations 
considered  in  [11]. 

We  can  then  apply  all  the  inversion  techniques  described  in  [10 ]  ,  [11  ] 
to  reconstruct  k(z)  or  O  (z)  .  However,  before  doing  so,  it  is  useful  to 

interpret  the  relations  (4.2).  We  note  first  that  if  <j>(z,r)  is  the  potential 

V 

of  the  earth  at  depth  z,  then  (J)(z,r)  is  related  to  <f>(z,t)  by 


f00 


4>(z,r)  = 


0 


<i>(z,X) J  (Ar)XdX 


r  v  dt 

<J>(z,t) 

J0 


.2  21/2  ' 
(t  +  r  ) 


(4.5) 


where  we  have  used  the  identity 


exp  -At  J  (Ar)dA  -  — - — 

0  (t  +  r  ) 


(4.6) 


The  transformation  (4.5)  has  an  interesting  property:  it  maps  solutions  of 


the  potential  equation  (2.6)  into  solutions  of 


s  y 

St 


d  S  V 

.z,t)  +  ~  £n0  (z)  v—  <f>(z,t)  =  0 

az  dz 


(4.7) 
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2  2  -1/2 

which  is  an  hyperbolic  equation.  To  see  this,  note  that  G^(r,t)  =  (t  +r  ) 


is  the  Green's  function  of  the  Laplacian,  i.e. 

Ag_  (t,r)  =  -  26  (t)  . 

0  r 

Then,  the  operator 

~A/8213  32  \  d  g  .  .  3 

A  =  I  — —  +  —  -5—  +  — —  j  +  q—  £na  (z)  ~— 
l  „  2  r  dr  ~  2  /  dz  dz 

\dz  dr  ' 


(4.8) 


(4.9) 


can  be  applied  to  both  sides  of  (4.5) ,  and  by  using  the  identity  (4.8) 
and  integrating  by  parts,  we  obtain 

r  .. 


A<j>(z,:r)  = 


0 


(W(f>  (z,  t)  )  Gq  ( t , r )  dt 


-  2(j)  (z,0) 


<5(r) 


3  V  x 

^4.(z,o)  - 


(4.10) 


where  W  denotes  the  perturbed  wave  operator 


W 


A  _  92  X 

\9z2  9t2/ 


2 1  lrCiz)  h 


(4.11) 


.3z  9t 

The  identity  (4.10)  shows  that  the  solutions  of  A<j)(z,r)  =  0  are  mapped  into 


solutions  of 

-V 

W(J>  (z,t)  =  0  , 


(4.12) 


with  the  initial  conditions 

(4.13) 

which  correspond  to  the  fact  that  the  system  (4.12)  is  originally  at  rest. 


V  3  V 

4>(z,0)  =  0,  T^j7<Hz,0)  =  0 


Fast  Cholesky  Recursions 


Then,  we  can  use  either  the  fast  Cholesky  or  Levinson  recursions  to 
solve  the  inverse  scattering  problem  associated  to  the  system  (4.4) .  However, 
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as  shown  in  Ill],  115],  in  order  to  apply  the  fast  Cholesky  recursions,  the 

V  V 

probing  waves  D(0,t)  and  U(0,t)  on  the  surface  of  the  medium  must  have  a 


very  specific  form,  i.e. 

>D(0,t)  =  6(.t)  +  d  (.0 ,  t)  1  (t) 

(4.14a) 

V 

U(0,t)  =  u(0, t) 1 (t) 

(4.14b) 

where  l(.t)  is  the  unit  step  function,  and  where  d(0,t)  and  u(0,t)  are 

V 

smooth  functions.  Thus,  D(0,t)  must  contain  a  leading  impulse  which  acts 


as  a  tag  indicating  the  wavefront  of  the  probing  wave. 

For  the  problem  considered  here,  the  potential  inside  the  earth 


can  be  expressed  as 


<Kz,r)  = 


2m  CO) 


.  ,  2  ,  2\ 1/2 
<(zz  +  r).  ' 


+  f(z, 


(4.15) 


0(0, X)  =  - - -  h(A)  .  (4.17b) 

2-na  '  (0) 

V 

Consequently,  if  Ah  (A)  =  Llh(.t)],  we  can  write 


V 

I 

V 

D(0,t) 

2tj01/2  (0) 

(6(t)  +  h(t)) 

(4.18a) 

V 

I 

V 

u(0,t) 

2itg1/2(0) 

h(t) 

(4.18b) 
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1/2  V  V 

so  that  modulo  multiplication  by  I/27TC  (0)  ,  the  waves  D(0,t)  and  U(0,t) 

V 

have  the  form  (4.14).  The  relation  d(0,t)  =  u(0,t)  =  h(t)  for  these 
waves  indicates  also  that  the  earth^s  surface  can  be  modeled  as  a  perfect 
reflector.  This  corresponds  to  the  fact  that  the  air  above  the  surface  of 
the  earth  acts  like  a  perfect  insulator. 

Then,  a  consequence  of  the  special  form  (4.14)  of  the  probing  waves 
is  that  the  waves  inside  the  scattering  medium  described  by  (4.4)  must 

have  the  form 

V 

D (z , t)  =  6(t-z)  +  d(z,t) 1 (t-z)  (4.19a) 

V 

U (z,t)  =  u(z,t) l(t-z)  .  (4.19b) 

By  substituting  (4.19)  inside  (4.4),  and  identifying  Coefficients  of  the 
impulse  6 (t-z)  on  both  sides  of  (4.4),  we  find  that 


fe  +  s)a<z-t) 
(fc-  ib)u(z't> 


k  (z)  u  (z ,  t) 


k(z)d(z,t) 


(4.20a) 


(4.20b) 


with 

k (z)  =  2u(z,  z+)  .  (4.21) 

After  discretization,  the  recursions  (4.20)  -  (4.21)  constitute  the 

fast  Cholesky  recursions  [11],  [15].  The  initial  data  for  these  recursions 

V 

is  d(0,t)  =  u(0,t)  =  h(t) .  The  relations  (4.20)  -  (4.21)  can  be  viewed  as 
using  a  layer- stripping  principle  to  identify  the  parameters  of  the  scattering 
medium.  Thus,  assume  that  the  waves  d(z,t)  and  u(z,t)  at  depth  z  have 
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been  computed.  The  reflectivity  function  k  (z I  is  obtained  from  (.4.21) 
and  is  used  in  (4.20)  to  compute  the  waves  d(z+A,  t)  and  u(z+A,  t)  at 
depth  z+A,  as  shown  in  Fig.  3.  The  effect  of  the  recursions  (4.20)  - 
(4.21)  is  therefore  to  identify  and  strip  away  the  layer  [z,  z+A).  Note 
that  the  Schur  recursions  of  Section  3  operated  according  to  a  similar 
principle. 

The  main  feature  of  the  fast  Cholesky  recursions  is  that  they  are  quite 
efficient:  let  L  be  the  maximum  depth  over  which  we  want  to  reconstruct 

the  medium,  and  let  A  =  L/N  be  the  step-size  which  is  used  to  discretize 

V 

the  fast  Cholesky  recursions.  Then,  by  observing  that  h(t)  needs  only  to 

be  known,  for  0  <_  t  <_  2L,  where  2L  is  the  two-way  travel  time  to  depth  L, 

and  computing  d(z,t)  and  u(z,t)  at  depth  z  only  for  0  <_  t  <_  2L  -  z,  we 

2 

find  Ill]  that  only  0 (N  )  operations  are  required  to  recover  k(z)  for 
0  <_  z_<  L.  In  addition,  it  was  shown  in  [24]  that  this  algorithm  is 
numerically  stable. 


Levinson  Recursions 


An  alternate  approach  is  to  formulate  the  inverse  scattering  problem  in 
terms  of  integral  equations.  Consider  the  Marchenko  integral  equations 


m11(z,t)  + 


ft  V 


h(t-x)m (z,T)dx  + 


fz  V 


h(t+T)m21(z,T)dx  =  0  (4.22a) 


V  ft  V 

h(z+t)  +  m21(z,t)  +  h(t+T)m11(z,x)dx 


ft  V 


+ 


h(t-x)m21(z,x)dx  =  0 


(4.22b) 
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with  -z  <  t  <  z.  Then,  it  is  shown  in  [11]  that  the  reflectivity  function 
k(z)  is  given  by 


k  (z)  =  - 

2m21(z,  z-)  , 

(4.23) 

and  that  m^Cz,*) 

and  m2^(z,»)  can  be  propagated  for 

increasing  values  of 

X 

z  by  using  the  Levinson  recursions 

/_3_  _3_\ 

\3z  +  3t / 

mH  (z , t)  =  ~  k(z)m21(z,t) 

(4.24a) 

(---) 
\9z  3t  / 

m21  (z,t)  =  “  k (z)m^1 (z, t) 

(4.24b) 

which  are  obtained  by  exploiting  the  Toeplitz  and  Hankel  structure  of  the 
kernels  appearing  in  (4.22).  The  initial  conditions  for  these  equations 
are 


mn  (O'  0)  =  m21(o,o)  =  0 


(4.25) 


and  in  the  propagation  of  (4,24)  we  use  the  boundary  conditions  m^(z,-z)  =  0 
and  (4.23),  where 


V  (z  V 

m21(z,z-)  =  -  h(2z)  -  hCz+T)m  (z,x)dx 

J_z 

•  z  V 

h(z-T)m21(z,T)dT  .  (4.26) 

J-z 

After  discretization,  the  Levinson  recursions  can  be  propagated  as  shown  in 

Fig.  4.  The  complexity  of  these  recursions  is  identical  to  that  of  the 

2 

fast  Cholesky  equations,  i.e.  they  require  0 (N  )  operations  to  reconstruct 
k(z)  for  0  <_  z  <_  L,  where  N  is  the  number  of  subintervals  which  are  used  to 
discretize  the  interval  [0,L] .  An  interesting  property  of  the  fast  Cholesky 
and  Levinson  recursions  (4.20)  and  (4.24)  is  that  they  have  the  same  form. 
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However,  the  support  of  the  functions  m.^(z,t)  and  m21(z,t)  is  -z  <_  t  <_  z, 

whereas  the  support  of  the  waves  d(z,t)  and  u(z,t)  is  z  t.  In  some  sense, 

the  Levinson  recursions  can  Be  viewed  as  being  the  complement  of  the  fast 

Cholesky  recursions:  they  rely  on  a  layer  accumulation  principle  where  at 

depth  z  we  identify  a  new  layer  and  accumulate  it  to  the  part  [0,z]  of  the 

medium  which  has  already  been  identified,  whereas  at  each  step  the  Schur 

and  fast  Cholesky  recursions  identify  and  strip  away  the  same  layer  from 

the  part  [z,°°)  of  the  medium  which  is  yet  to  be  identified.  An  additional 

difference  between  the  fast  Cholesky  and  Levinson  recursions  is  that  the 

fast  Cholesky  recursions  correspond  to  an  initial  value  problem  where 

all  the  information  about  the  medium  is  contained  in  the  initial  conditions 

d(0,t)  and  u(0,t),  while  for  the  Levinson  recursions  the  identification 

of  a  new  layer  requires  at  every  step  the  evaluation  of  the  integral  (4.26) 

V 

where  the  xnformatxon  about  the  medium  is  contained  in  h (t ) . 

It  turns  out  that  the  above  reconstruction  procedure  for  the  earth's 
resistivity  is  not  new,  and  appears  in  disguised  form  in  Kunetz  and 
Rocroi  .  [  8]  for  the  case  of  a  discrete  medium  with  layers  of  equal  thickness. 
However,  Kunetz  and  Rocroi  did  not  identify  the  recursions  that  they  obtained 
as  the  Levinson  recursions. 

The  previous  reconstruction  procedure  can  also  be  related  to  that  of 
Coen  and  Yu  [6]  by  noting  that  if 

A (z , t)  =mil(z,t)  +  m21(z,t)  ,  (4.27) 

then  A(z , t)  satisfies  the  integral  equation  (see  [11]) 

V  fz  V  v 

h (z+t)  +  A (z , t)  +  (h(t-T)  +  h(t+T) ) A(z,T)dT  =  0  (4.28) 

■'-z 
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with  -z  <_  t  <_  z,  which  is  the  equation  used  by  Coen  and  Yu.  Then 

1/2  (Z 

=  1  +  A (z , t) dt  (4.29) 

•'-z 

so  that  a(z)  can  be  reconstructed  directly  from  A(*,-).  The  advantage  of 
this  method  over  the  procedure  described  above  is  that  since  0  ( . )  is 
smoother  than  k(-),  it  is  easier  to  reconstruct.  However,  Coen  and  Yu 
were  unaware  of  the  existence  of  a  fast  algorithm  to  solve  the  integral 
equation  (4.28). 
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5.  Interpretation  and  Computation  of  the  Inversion  Data 


The  inversion  procedures  described  above  rely  on  R  (A)  ,  or  equivalently 

Lj 

V 

on  Slichter's  kernel  K(A)  ,  and  on  the  function  h (t)  to  reconstruct  the  earth's 
resistivity.  But  the  given  data  is  the  apparent  resistivity  p  (r)  obtained 

3 

from  the  Schlumberger  electrode  configuration.  The  problem  of  computing  K(A) 

from  p  (r)  was  solved  by  Ghosh  [25]  who  used  the  expression 

3. 


K(A)  =  0(0) 


f°° 

—  p  (r)  J  (Ar)  dr 
r  a  1 

J0 


(5.1) 


which  is  obtained  by  combining  (2.16)  and  (2.30),  so  that 


K(A) 


2tto  (0) 
I 


A4>  (o,A) 


(5.2) 


and  by  using  the  identity  (2.10)  of  Hankel  transforms  and  the  definition 
(2.1)  of  the  apparent  resistivity.  Then,  by  substituting 


x 

r  =  e  , 


A  = 


(5.3) 


inside  (5.1)  and  denoting 

p  (x)  =  p  (eX) ,  K(y)  =  K(e~Y)  , 

3  3 


we  obtain  the  convolution  integral 


K(y) 


r  ~ 

p  (x)j  (exp(- (y-x) ) )dx 

J0  a  1 


(5.4) 


(5.5) 


which  can  be  implemented  by  discrete  convolution  techniques  [26] . 
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V 

The  problem  of  computing  h(t)  from  p  (r)  is  more  difficult.  The 

a 

V 

first  step  is  to  obtain  a  physical  interpretation  of  h(t).  From  (4.15),  we 
find  that 


<KO,r) 


I 

2ttc(0) 


dt 


.  2  2  1/2 
(t  +  r  )  ' 


(5.6) 


The  first  term  in  this  expression  is  the  potential  associated  to  a  homogeneous 

earth  with  conductivity  a(0),  and  the  second  term  describes  the  effect  of 

inhomogeneities  in  the  earth's  resistivity.  However,  to  describe  the  potential 

on  the  surface  of  the  earth,  instead  of  assuming  that  the  earth's  resistivity 

is  inhomogeneous  and  that  a  single  current  source  is  located  at  the  origin 

of  coordinates,  by  Maxwell's  method  of  images  ([27],  [5],  p.  197)  we  can 

assume  that  the  earth  is  homogeneous  with  conductivity  0 (0) ,  but  that  some 

additional  fictitious  current  sources  have  been  added  on  the  vertical  axis . 

V 

In  this  case,  if  h(t)dt  is  the  strength,  relative  to  the  strength  I  of  the 
actual  current  source,  of  a  source  located  at  depth  t  along  an  infinitesimal 
segment  of  length  dt,  the  potential  created  at  the  point  (0,r)  on  the  surface 
of  the  earth  is 


V 

I  h  (t)  dt 

2TO  (0)  ,  2  2  1/2 

(t  +  r  ) 


(5.7) 


Note  that  in  order  to  guarantee  that  the  vertical  component  of  the  current 

density  created  by  the  fictitious  sources  is  zero  on  the  surface  of  the  earth, 

V 

the  function  h(t)  must  be  symmetric  with  respect  to  the  origin,  i.e.  sources 


must  be  located  above  the  surface  of  the  earth  as  well  as  below.  By  superposition, 
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V 

the  function  h(t)  appearing  in  (5.6)  can  therefore  be  viewed  as  the  fictitious 
current  source  profile  equivalent  to  the  inhomogeneous  conductivity  profile 
0(z)  . 

V 

The  function  h(t)  was  the  starting  point  of  the  inversion  method  of 
Kunetz  and  Rocroi  [ 8 ] .  However,  it  is  not  as  easy  to  compute  this  function 
from  the  potential  cf>(0,r)  or  the  apparent  resistivity  p  (r)  as  it  appears, 

cl 

V 

To  see  why  this  is  so,  note  from  (4.16)  that  in  order  to  obtain  h(t)  from 
<j)(0,r),  we  need  first  to  compute  the  Hankel  transform  4>(0,A)  followed  by  an 
inverse  Laplace  transform.  But  inverse  Laplace  transforms  are  hard  to 
implement.  Instead,  it  is  preferable  to  discretize  the  integral  equation  (5.6) 
and  to  solve  the  resulting  system  of  linear  equations.  In  terms  of  p  (r) ,  we 

cL 

find  from  (2.1)  that 


P  (r) 

3. 


( 

0(0)  \ 


1  +  2r~ 


3  r  y 


h  (t) 


0 


dt  ^ 

(t2  +  rV/2 


(5.8) 


(see  18],  [5],  Chapter  10),  which  can  also  be  discretized  and  inverted. 

An  alternate  method  of  computing  Y(t) ,  which  was  proposed  by  Kunetz  and 

V 

Rocroi  [8  ],  is  to  denote  by  H(k)  the  Fourier  transform  of  h(t)  and  to  introduce 
the  spectral  density  function 

W(k)  =  1  +  H (k)  +  H(-k)  .  (5.9) 


Then,  by  observing  from  (4.18)  that  the  left  reflection  coefficient  of  the 

two-component  scattering  system  is 

R_  (jk)  =  H(k)/(1  +  H (k) )  (5.10) 

L 

and  is  bounded  by  one,  i.e.  ) R  ( jk) |  <_  1,  we  can  conclude  that  the  spectral 

Li 

density  W(k)  is  positive,  i.e.  W(k)  0  for  all  k.  By  using  the  integral 
equation  (5.8),  we  also  find  that 
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P  (r) 

ci 


2r' 


a(o)iT 


,CO 

Jo 


W(k)K  (kr)kdk 


(5.11) 


where  K  (.*)  is  the  modified  Bessel  function  of  order  one.  The  problem  of 

inverting  the  modified  Hankel  transform  (.5.11)  is  analogous  to  that  of  inverting 

a  Laplace  transform,  but  by  using  the  positivity  of  W(k) ,  Kunetz  and  Rocroi 

were  able  to  formulate  the  inversion  of  (5.11)  for  a  discrete  set  of  sampled 

values  of  r  as  a  quadratic  programming  problem.  Then,  given  the  reconstructed 
V  I  I 

W(k) ,  h ( I t [ )  is  the  inverse  Fourier  transform  of  W(k)  -  1. 


-31- 


6.  Conclusion 

In  this  paper,  we  have  considered  the  problem  of  reconstructing  the 
resistivity  profile  of  a  layered  earth  probed  by  direct  current  from 
potential  measurements  on  the  surface  of  the  earth.  It  was  shown  that  this 
problem  could  be  transformed  into  an  equivalent  inverse  scattering  problem, 
to  which  efficient  signal  processing  algorithms  such  as  the  Schur,  fast 
Cholesky  and  Levinson  recursions  can  be  applied.  These  algorithms  recon¬ 
struct  the  resistivity  of  the  earth  layer  by  layer  in  a  recursive  way,  and 
require  only  a  small  number  of  operations.  In  this  context,  it  was  shown 
that  the  recursions  obtained  by  Pekeris  13]  and  Koefoed  [4] ,  [5]  for 
recovering  the  resistivity  of  the  earth  were  identical  to  the  discrete 
Schur  recursions,  and  that  the  reconstruction  method  of  Kunetz  and 
Rocroi  18]  was  actually  based  on  the  Levinson  recursions. 

One  difficulty  associated  with  these  reconstruction  methods  is  that 
they  do  not  operate  directly  on  the  given  data,  which  is  the  apparent 

resistivity  of  the  earth,  but  on  Slichter's  kernel  K(X) ,  or  on  the 

V 

fictitious  current  source  profile  h(t).  equivalent  to  the  inhomogeneous 
conductivity  profile  a( z) .  Efficient  convolution  techniques  exist  to 

compute  Slichter's  kernel  from  the  apparent  resistivity,  but  the  problem 

V  V 

of  computing  h(t)  is  more  difficult.  No  efficient  method  of  obtaining  h(t) 

exists,  short  of  brute  force  discretization  of  the  integral  equation  satisfied 

V 

by  h(t).  This  problem  deserves  therefore  further  attention.  Another 
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topic  of  research,  which  is  currently  under  investigation,  is  the  study  of 
the  numerical  behavior  of  the  algorithms  described  above  when  they  operate 
on  synthetic  or  real  data.  The  fast  Cholesky  and  Levinson  recursions 

are  known  to  be  stable,  but  the  addition  of  noise,  or  imperfections  in 
the  data  due  to  bandlimitations ,  can  degrade  the  performance  of  these 
algorithms.  It  would  therefore  be  desirable  to  develop  inversion  techniques 
which  can  incorporate  a  priori  information  on  the  resistivity  profile  and 
on  the  noise  level. 
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FIGURE  CAPTIONS 


Fig.  1 
Fig.  2 

Fig.  3 

Fig.  4 


:  Elementary  filter  sections  associated  to  the  two-component  System. 

:  The  aggregate  medium  obtained  by  composing  the  elementary  filter 

sections . 

:  (a)  Propagation  of  d(z,t)j  and  b)  propagation  of  u(z,t)  via  the 

fast  Cholesky  recursions. 

:  (a)  Propagation  of  m^Czjt);  and  b)  propagation  of  m2^(z,t)  with 

the  Levinson  recursions. 
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